TP 3 - Partie 2

Etudiants : Doré Maxime, Pillard Cécile, Bennis Zineb, Debbagh Boutarbouch Yasmine, Jacquiet Morgane

defaultW <- getOption("warn")
options(warn = -1)


library(lhs)
library(sensitivity)
library(plotly)
library(DiceKriging)
library(randtoolbox)
library(DiceDesign)

options(warn = defaultW)

1. Calcul analytique des indices de Sobol

#Cette fonction prend en argument une matrice X qui contient 3 colonnes de longueur n
# des valeurs réparties uniformément entre -pi et pi
#Les colonnes correspondent respectivement à x1,x2 et x3
ishigami <- function(X, coeff){
    return(sin(X[,1]) + coeff[1]*sin(X[,2])^2 + coeff[2]*X[,3]^4*sin(X[,1]))
}

#Cette fonction calcule les indices de Sobol Analytiques 
# du premier et du second ordre pour la fonction Ishigami
ishigami.sobol.indices <- function(coeff){
    
    V = 1/2 + coeff[1]^2/8 + ((coeff[2]^2)*pi^8)/18 + (coeff[2]*pi^4)/5
    V1 = (1/2)*(1+coeff[2]*((pi^4)/5))^2
    V2 = (coeff[1]^2)/8
    V3 = 0
    V12 = 0
    V13 = coeff[2]^2*(pi^8)*(8/225) 
    V23 = 0
    
    # Initialisation des coeffs
    S1 <- St <- array(0,3)
    
    # Coeffs premier ordre
    S1[1] = V1/V
    S1[2] = V2/V
    S1[3] = V3/V #=0
    
    # Coeffs second ordre
    S12 = V12/V #=0
    S13 = V13/V
    S23 = V23/V #=0
    S2 <- c(S12, S23, S13)
    
    # Coeffs totaux
    St[1] = (V1+V13)/V
    St[2] = V2/V
    St[3] = V13/V
    list(Coeff_Ordre1=S1, Coeff_Ordre2=S2, Coeff_Ordre_Tot=St)
}

n = 4000

X = matrix(NA, n, 3) #Creation d'une matrice vide

x1 = runif(n,-pi,pi)
x2 = runif(n,-pi,pi)
x3 = runif(n,-pi,pi)

X[,1] = x1
X[,2] = x2
X[,3] = x3

#initialisation des coefficients A et B de la fonciton d'Ishigami:
coeff = c(7, 0.1)

cat("Indices de Sobol pour la fonction Ishigami : \n")
(ishigami.sobol.indices(coeff))

Indices de Sobol pour la fonction Ishigami : 

$Coeff_Ordre1

1. 0.313905191147811
2. 0.442411144790041
3. 0

$Coeff_Ordre2

1. 0
2. 0
3. 0.243683664062148

$Coeff_Ordre_Tot

1. 0.557588855209959
2. 0.442411144790041
3. 0.243683664062148

2. Simulation des effets principaux et effets d'interaction

# fonctions théoriques 
g0 = coeff[1]/2

g1 = sin(sort(x1))*(1 + coeff[2]/5 * pi**4)
g2 = coeff[1] * (sin(sort(x2))**2 - 1/2)
g3 = rep(0, length(x3))

g12 = rep(0, length(x1))
g13 = coeff[2] * sin(sort(x1)) * (sort(x3)**4 - 1/5 * pi**4)
g23 = rep(0, length(x3))

Visualisation des interactions du premier ordre

ish = ishigami(X, coeff)
# On centre ish
ish_centre = ish - g0

# polynomes locaux
g1_loc = predict(loess(ish_centre ~ x1, span=0.5, family="gaussian"))
g2_loc = predict(loess(ish_centre ~ x2, span=0.5, family="gaussian"))
g3_loc = predict(loess(ish_centre ~ x3, span=0.5, family="gaussian"))

options(repr.plot.width = 14, repr.plot.height = 5) 
par(mfrow=c(1,3))

plot(x1, ish_centre, type = 'p')
# Prediction
lines(sort(x1), g1_loc[order(x1)], lwd=6, col="blue")
# Theorique
lines(sort(x1), g1, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)

plot(x2, ish_centre, type = 'p')
lines(sort(x2), g2_loc[order(x2)], lwd=6, col="blue")
lines(sort(x2), g2, lwd=3, col="red")
legend(x='top', legend=c("Estimé", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)

plot(x3, ish_centre, type = 'p')
lines(sort(x3), g3_loc[order(x3)], lwd=6, col="blue")
lines(sort(x3), g3, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)

Visualisation des interactions du second ordre

On se rend compte qu'on ne peut pas si facilement reproduire la méthode précédente dans les cas d'intéractions du second ordre. On va donc utiliser une méthode un peu plus "manuelle". Par exemple pour simuler g13 :

• On fixe X1
• On fixe X3
• On regarde l'évolution de Ishigami selon X2

Vu que g13 est la seule fonction non nulle, on ne représentera que celle-ci.

g13_chap = c()
x1 = c()
x3 = c()
nb_point = 300
#Initialisation des valeurs fixées prédéfinies
X1_grid <- X3_grid <- seq(from = -pi,to=pi,length.out =nb_point)

#Voir avancement du programme
c=0 
pb <- txtProgressBar(min = 1, max = length(X1_grid), style = 3)

for (i in X1_grid) {
    c=c+1
    Sys.sleep(0.001)
    setTxtProgressBar(pb, c)
    #on fixe X1
    X1 = rep(i,nb_point)
    for (j in X3_grid) {
        #On fixe X3 
        X3 = rep(j,nb_point)
        #on définit X2 qu'on fait varier
        X2 = runif(nb_point, min=-pi, max=pi)
        contenu = c(X1,X2,X3)
        #on remplit la matrice pour X1 et X3 fixés
        X_fixe = matrix(contenu, nrow = nb_point, ncol=3)
        Ish_X_fixe = ishigami(X_fixe, coeff)
        g1_aux = sin(sort(X1))*(1 + coeff[2]/5 * pi**4)
        g2_aux = coeff[1] * (sin(sort(X2))**2 - 1/2)
        g13_chap_ij = mean(Ish_X_fixe)-g0-g1_aux[1]-g2_aux[1]
        g13_chap = c(g13_chap, g13_chap_ij)  
        #pour l'affichage
        x1 = c(x1, i)
        x3 = c(x3, j)
    }
}

  |======================================================================| 100%

defaultW <- getOption("warn")
options(warn = -1)

fig <- plot_ly(x = ~X[,1], y = ~X[,3], z = ~ish_centre, mode='markers', marker = list(symbol = 'circle', sizemode = 'diameter', size = 2), type='scatter3d') 
fig <- fig %>% add_trace(x = ~x1, 
                         y = ~x3, 
                         z = ~g13_chap, type = 'mesh3d')
fig

options(warn = defaultW)

Si l'image ci-dessus ne s'affiche pas, cf l'image .png envoyé dans le mail

3. Estimation des indices de Sobol par métamodèle processus gaussien (krigeage)

Construction du plan d'expérience

Quel type de plan d’expériences utiliseriez-vous et pourquoi ?

Les propriétés attendues d'un plan d'expériences numériques sont :

• la couverture de l'espace des paramètres d'entrées (=discrepance faible)
• la robustesse en sous projection
• la séquentialité.

Nous allons calculer la discrepance pour différents plan d'expériences.

critere=0 
for (i in 1:10000)
{
  x = matrix(runif(n=2*9),nrow=9,ncol=2) # Plan pur MonteCarlo
  if (mindist(x) > critere) # selection du plan courant qui maximise 
                            # la distance minimale entre 2 de ses points
  { 
    x1 = x  # Plan Maximin
    critere = mindist(x)
  }
}


dis_MC = discrepancyCriteria(x,type=c('M2','C2')) 
cat("Critères discrépance M2 et C2 pour le plan pur Monte Carlo = ", 
    dis_MC$DisM2, " - ", dis_MC$DisC2,"\n")



dis_Maximin = discrepancyCriteria(x1,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour le plan Maximin = ", 
    dis_Maximin$DisM2, " - ", dis_Maximin$DisC2,"\n")


x2 = factDesign(3,3)
dis_Fact = discrepancyCriteria(x2$design,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour le plan Factoriel = ", 
    dis_Fact $DisM2, " - ", dis_Fact $DisC2,"\n")


# Création d'une suite de Sobol de 9 points en dimension 3
x3 = sobol(9,3)
dis_Sobol = discrepancyCriteria(x3,type=c('M2','C2'))
cat("Critères discrépance M2 et C2 pour la suite de Sobol = ", 
    dis_Sobol$DisM2, " - ", dis_Sobol$DisC2,"\n")

Critères discrépance M2 et C2 pour le plan pur Monte Carlo =  0.2656244  -  0.2316414 
Critères discrépance M2 et C2 pour le plan Maximin =  0.1159467  -  0.1103984 
Critères discrépance M2 et C2 pour le plan Factoriel =  0.3828553  -  0.3167492 
Critères discrépance M2 et C2 pour la suite de Sobol =  0.1644181  -  0.1305912 

Le plan à discrepance faible avec une suite de Sobol nous permet de minimiser le critère de discrepance. De plus, ce plan permet la séquentialité car il n'est pas nécessaire de recalculer le plan total lorsque l'on rajoute des nouveaux points. C'est donc ce plan que nous utiliserions.

Construire cette base d’apprentissage à partir d’un un plan LHS (optimisé) de taille N = 100 points

# Nombre d'évaluations de f
N = 100

pex = (optimumLHS(N, 3)-0.5)*2*pi
options(repr.plot.width = 5, repr.plot.height = 5) 
plot(pex, main="LHS")

# Evaluation f selon le plan pex
pex <- data.frame(pex)
y <- ishigami(pex, coeff)

Estimer un métamodèle PG et contrôler sa qualité de prédiction

# Modèle de Krigeage

m.krigeage <- km(design=pex, response=y,control=list(trace=FALSE))
options(repr.plot.width = 5, repr.plot.height = 10) 
plot(m.krigeage)

Evaluation de sa qualité de prédiction

#Définition de la fonction qui calcule le Q2 : 
q2 <- function(y, y_pred){
    return ( 1- ( (sum((y-y_pred)^2)) / sum((mean(y)-y_pred)^2) ) )
}

cross_validation <- function(pex,nb_fold){
    # Evaluation de la capacité de prédiction par validation croisée

    #Randomly shuffle the data
    pex_sample<-pex[sample(nrow(pex)),]  
    #Create 10 equally size folds
    folds <- cut(seq(1,nrow(pex_sample)),breaks=nb_fold,labels=FALSE)

    #Perform k fold cross validation
    xval_err = 0
    Q2 = 0
    vec_y_test = c()
    vec_y_test_pred = c()
    for(i in 1:nb_fold){
        #Séparation des données pour la cross validation
        testIndexes <- which(folds==i,arr.ind=TRUE)
        testData <- pex_sample[testIndexes, ]
        trainData <- pex_sample[-testIndexes, ]
        
        #Calcul du modèle sur le jeu train
        y_train = ishigami(trainData, coeff)
        m.krigeage <- km(design=trainData, response=y_train,control=list(trace=FALSE))
        
        #prediction sur le jeu de test
        y_test_pred = predict.km(m.krigeage, testData, "UK", se.compute=FALSE, 
                                 checkNames=FALSE)$mean
        vec_y_test_pred = c(vec_y_test_pred,y_test_pred)

        y_test = ishigami(testData, coeff)
        vec_y_test = c(vec_y_test,y_test)


        e = y_test_pred-y_test
        xval_err = xval_err + e%*%e/length(y_test) 

        Q2 = Q2 + q2(y_test, y_test_pred)
    }

    rmse = sqrt(xval_err/nb_fold)
    Q2 = Q2/nb_fold  
    
    return(list(RMSE = rmse, Q2=Q2, Vec_pred = vec_y_test_pred, Vec_reel = vec_y_test))
    
}

list_result = cross_validation(pex,10)

cat("Valeur de la RMSE : " , list_result$RMSE) #doit être le plus petit possible

Valeur de la RMSE :  1.304057

cat("Valeur du Q2 : " ,list_result$Q2) # doit etre le plus proche  de 1

Valeur du Q2 :  0.8491658

#Graphe des valeurs prédites VS valeurs réelles :
options(repr.plot.width = 5, repr.plot.height = 5) 

y_test = list_result$Vec_reel
y_test_pred = list_result$Vec_pred
plot(y_test, y_test_pred, main = "Valeurs prédites vs Valeurs réelles", xlab="Reelles", ylab="Predites")
reg <- lm(y_test_pred ~ y_test)
abline(reg, col='red')
lines(y_test,y_test, col='blue')
legend(x='bottom', legend=c("Reg Predite", "Reg Reelle"), col = c('red','blue'), cex= 0.9, lty=1, bty='n')

Estimer les indices de sobol du 1er ordre et totaux

m.krigeage <- km(design=pex, response=y,control=list(trace=FALSE))

kriging.mean <- function(Xnew, m){
  predict.km(m, Xnew, "UK", se.compute = FALSE, checkNames = FALSE, control=list(trace=FALSE))$mean
}


SA.metamodel <- fast99(model = kriging.mean, factors = 3, n = 1000,
                       q = "qunif", q.arg = list(min = -pi, max = pi), m = m.krigeage)


#SA_test <- fast99(model = ishigami.fun, factors = 3, n = 1000,
           # q = "qunif", q.arg = list(min = -pi, max = pi))

options(repr.plot.width = 5, repr.plot.height = 5) 
plot(SA.metamodel, cex=0.8)

Interprétation :

La variable qui a le plus d’influence sur la variance de la sortie (c'est à dire avec l'indice total le plus élevé) est la variable X1, avec un indice total d'environ 60% et un indice du premier ordre d'environ 30%. (X1 explique à elle seule environ 30% de la variance de Y)

La variable X2 n’intervient que seule dans la fonction Ishigami, cela se retrouve sur le graphique ci-dessus où son indice du premier ordre est quasiment équivalent à son indice total. (environ 40% de la variance de Y est expliquée par X2)

La variable X3 n’a aucune influence seule, mais a une influence relativement importante en intéraction et notamment avec X1 (comme on le retrouve d'ailleurs dans la formule d'Ishigami).

On en déduit que la variance de Y est due pour 40% à X2, 30% à X1 et 30% à l’interaction entre X1 et X3.

cat("Indices de Sobol pour la fonction Ishigami : \n")
(ishigami.sobol.indices(coeff))
print(SA.metamodel)

Indices de Sobol pour la fonction Ishigami : 

$Coeff_Ordre1

1. 0.313905191147811
2. 0.442411144790041
3. 0

$Coeff_Ordre2

1. 0
2. 0
3. 0.243683664062148

$Coeff_Ordre_Tot

1. 0.557588855209959
2. 0.442411144790041
3. 0.243683664062148

Call:
fast99(model = kriging.mean, factors = 3, n = 1000, q = "qunif",     q.arg = list(min = -pi, max = pi), m = m.krigeage)

Model runs: 3000 

Estimations of the indices:
   first order total order
X1 0.296888140   0.5186778
X2 0.455851334   0.5031961
X3 0.004771931   0.2350741

coeff_ordre1_metam=SA.metamodel$D1/SA.metamodel$V
coeff_tot_metam = 1-((SA.metamodel$Dt)/SA.metamodel$V)

coeff_ordre1=ishigami.sobol.indices(coeff)$Coeff_Ordre1
coeff_ordre_tot=ishigami.sobol.indices(coeff)$Coeff_Ordre_Tot



cat("Différence absolue les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre 1 :\n")
cat(abs(coeff_ordre1_metam-coeff_ordre1), '\n \n')
cat("Différence les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre total :\n")
cat(abs(coeff_tot_metam-coeff_ordre_tot))

Différence absolue les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre 1 :
0.01701705 0.01344019 0.004771931 
 
Différence les valeurs estimées avec le métamodèle et les valeurs théoriques des coefficients d'ordre total :
0.03891105 0.06078491 0.008609577

Estimer les effets principaux en utilisant le métamodèle par simulation intensive

n = 4000
X = matrix(NA, n, 3) #Creation d'une matrice vide

x1 = runif(n,-pi,pi)
x2 = runif(n,-pi,pi)
x3 = runif(n,-pi,pi)

X[,1] = x1
X[,2] = x2
X[,3] = x3

#initialisation des coefficients A et B de la fonciton d'Ishigami:
coeff = c(7, 0.1)

ish_krig = kriging.mean(X,m.krigeage)
# On centre ish_krig
ish_krig_centre = ish_krig - g0

# polynomes locaux
g1_loc_krig = predict(loess(ish_krig_centre ~ x1, span=0.5, family="gaussian"))
g2_loc_krig = predict(loess(ish_krig_centre ~ x2, span=0.5, family="gaussian"))
g3_loc_krig = predict(loess(ish_krig_centre ~ x3, span=0.5, family="gaussian"))

options(repr.plot.width = 14, repr.plot.height = 5) 
par(mfrow=c(1,3))

plot(x1, ish_krig_centre, type = 'p')
# Prediction
lines(sort(x1), g1_loc_krig[order(x1)], lwd=6, col="blue")
# Theorique
lines(sort(x1), g1, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé (krigeage)", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)

plot(x2, ish_krig_centre, type = 'p')
lines(sort(x2), g2_loc_krig[order(x2)], lwd=6, col="blue")
lines(sort(x2), g2, lwd=3, col="red")
legend(x='top', legend=c("Estimé (krigeage)", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)

plot(x3, ish_krig_centre, type = 'p')
lines(sort(x3), g3_loc_krig[order(x3)], lwd=6, col="blue")
lines(sort(x3), g3, lwd=3, col="red")
legend(x='topleft', legend=c("Estimé (krigeage)", "Theorique"), col=c("blue","red"), lwd=c(2,1),lty=1,text.font=1.5,bty='n',cex=1.4)

Avec seulement 100 simulations, on a des bonnes approximations. Cela montre l'intérêt de passer par des métamodèles.

Refaire pour plusieurs tirages aléatoires différents

# Nombre d'évaluations de f
N_seq = c(20,50,100,150)
nb_eval = 5
nb_fold = 10

df_rmse <- list()
df_Q2 <- list()

i= 1
for (Ni in N_seq) {

    vec_rmse = c()
    vec_Q2 = c()
    for (k in 1:nb_eval){
        pex_i = (optimumLHS(Ni, 3)-0.5)*2*pi
        # Evaluation f selon le plan pex
        pex_i <- data.frame(pex_i)
        list_result_i = cross_validation(pex_i,10)
        vec_rmse <- c(vec_rmse, list_result_i$RMSE[1])
        vec_Q2 <- c(vec_Q2, list_result_i$Q2)  
    }
    print(vec_rmse)
    df_rmse[[i]] <- vec_rmse
    df_Q2[[i]] <- vec_Q2
    i = i+1
}

#options(repr.plot.width = 5, repr.plot.height = 5) 
#plot(pex, main="LHS")

[1] 4.594274 4.969337 3.978402 3.344505 3.276791
[1] 2.950537 3.295803 2.225401 2.274905 3.703246
[1] 1.375695 1.154038 1.461013 1.339915 1.335368
[1] 0.5825002 0.8992492 0.8225002 0.6376081 0.7362206

boxplot(df_rmse,main='Boxplot de RMSE selon la taille du LHS',
        names=c(20,50,100,150),
        xlab='Taille LHS',ylab='RMSE',col='lightblue')

Plus le nombre de simulations est élevé, plus la RMSE diminue.